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ABSTRACT 

For spectroscopically confirmed type la supernovae we evaluate models of intrinsic brightness vari- 
ations with detailed data/Monte-Carlo comparisons of the dispersion in the following quantities: 
Hubblc-diagram scatter, color difference (B — V — c) between the true B — V color and the fit- 
ted color (c) from the SALT-ii light curve model, and photometric rcdshift residual. The data sample 
includes 251 ugriz light curves from the 3-season Sloan Digital Sky Survey -II, and 191 griz light 
curves from the Supernova Legacy Survey 3-year data release. We find that the simplest model of 
a wavelength independent (coherent) scatter is not adequate, and that to describe the data the in- 
trinsic scatter model must have wavelength-dependent variations resulting in a ~ 0.02 mag scatter 
in B — V — c. Relatively weak constraints are obtained on the nature of intrinsic scatter because 
a variety of different models can reasonably describe this photometric data sample. We use Monte 
Carlo simulations to examine the standard approach of adding a coherent scatter term in quadrature 
to the distance-modulus uncertainty in order to bring the reduced x 2 to unity when fitting a Hubble 
diagram. If the light curve fits include model uncertainties with the correct wavelength dependence of 
the scatter, we find that this approach is valid and that the bias on the dark energy equation of state 
parameter w is much smaller (~ 0.001) than current systematic uncertainties. However, incorrect 
model uncertainties can lead to a significant bias on the distance moduli, with up to ~ 0.05 mag 
redshift-dependent variation. This bias is roughly r educed in half after applying a Malmquist bias 
correction. For the recent SNLS3 cosmology results (jSullivan et al.ll2~011f ) we estimate that this effect 
introduces an additional systematic uncertainty on to of ~ 0.02, well below the total uncertainty. 
However, this uncertainty depends on the samples used, and thus this small w-uncertainty is not 
guaranteed in future cosmology results. For example, the u;- uncertainty for SDSS+SNLS (dropping 
the nearby SNe) increases to ~ 0.04. 
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1. INTRODUCTION 

For more than a decade Type la supernovae (SN la) 
have been used as standard candles to measure lu- 
minosity distances. These distances, along with 
the associated rcdshifts, have be en used to mea- 
sure properties of d a rk energy dRiess et al.l 119981: 
Perlmutter et al.l 119991: iRiess et all 12004 lAstier et all 
20061: iWood-Vasev et all l2007t iFreedman et all 120091: 
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IKessler et all 1200911 IConlev et all l20TTh . The uncor- 
rected variation in the SN la peak brightness is ~ 1 
magnitude, and this variation is reduced to ~ 0.1 mag af- 
ter empiric al corrections based on the measured stretch 
and color (|Phillipsl Il99l . This 0.1 mag intrinsic scat- 
ter increases the scatter in the Hubble diagram well be- 
yond what is expected from the distance modulus uncer- 
tainties, and the resulting cosmology fits have reduced 
chi-squared (x 2 = X 2 /-^dof) significantly larger than 
unity. To obtain ^ = 1, all SNIa- cosmology analyses 
to date have introduced an ad-hoc intrinsic scatter term 
(oint ~ 0.1 mag) that is added in quadrature to the mea- 
sured distance moduli uncertainties. 

This procedure of adding a constant ad-hoc scatter 
term would be correct if the unknown source of in- 
trinsic variation is independent of redshift, and if it is 
fully coherent such that the variation is the same for all 
wavelengths and passbands. IKessler et all (|2010l here- 
after K10) found evidence contradicting a coherent vari- 
ation in a study comparing the photo- z precision be- 
tween data and Monte Carlo simulations (MC). Using 
the coherent-scatter model in the MC underestimated 
the fitted photo- z precision observed in the data, while 
simulating a model us i ng co lor variations gave better 
agreement. I Guv et all (|2010l hereafter G10) examined 
residuals from the SALT-ii training and showed that the 
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variation about the best-fit spectral model is wavelength 
dependent; this wavelength-dependent u ncertainty is in- 
cluded in the light curve fitting model. [Marrinc r et al.l 
(|201lL hereafter Mf f ) presented a more formal treatment 
of (Ti n t based on an intrinsic scatter covariance matrix 
that depends on a coherent term, a stretch term, and a 
color termJ3 They compared SDSS-II SN la data and 
simulations within this framework and suggest that the 
coherent and color terms are significant, while the stretch 
term is negligible. 

The importance of understanding the nature of intrin- 
sic scatter is tied to understanding systematic uncertain- 
ties in cosmology analyses using SNc la. If this scatter 
is truly random, as suggested by explosion models show- 
ing b rightness variation with viewing angle (|Kasen et al.l 
120091 hereafter KRW09), then there is no intrinsic bias 
and the uncertainty will decrease with increasing sam- 
ple size. Even in this optimistic scenario a wavelength 
dependent scatter results in a redshift dependent disper- 
sion that must be properly accounted for. If this scatter 
depends on more subtle physics related to the explosion 
mechanism and the host galaxy environment, there could 
be additional redshift-dependent effects not yet detected 
with current data samples, but that become apparent in 
future surveys with much larger samples. 

Significant effort to reduce this scatter has been at- 
tempted using near in frared (NIR) p h otom etry and spec- 
troscopic features. iMandel et al.1 (|2011| ) report that 
optical+NIR photometry result in a Hubble scatter 
that is ~ 30% smaller compared to using only opti- 
cal data. A decade of effort on spectroscopic correla- 
tions can be summar ized with results from three groups. 
iBlondin et al.l (|2011f ) examined spectra from the CfA 
Supernova Program and used spectral features to re- 
duce the Hubble scatter by at most 10%, albeit with 
only 2a significance. F rom the Berkeley SNIa Program, 
iSilverman et al.l (|2012t) examined 108 high quality SNe la 
with a spectrum taken within 5 days of maximum bright- 
ness and found similar result s. The best Hubbl e scat- 
ter reduction was obtained by iBailev et al.l ()2009l ) using 
very high quality sp ectra from the Supernova Factory 
(jAldering et al.ll2002| ). Using spectra within 2.5 days of 
peak brightness, they scanned every possible flux-ratio in 
~ 40 A bins and found a minimum Hubble scatter using 
the ratio F(642 nm)/F(443 nm); the resulting scatter is 
about 25% lower compared to the traditional photomet- 
ric analysis with the SALT-ii light curve model. 

To realize significant reductions (~ 30%) in the Hub- 
ble scatter requires optical photometry combined with 
either rest-frame NIR photometry or very high quality 
spectra near the epoch of peak brightness. Both of these 
supplemental data samples are difficult to obtain at low 
redshift, and it is not yet clear what resources could be 
allocated to obtain large data samples for higher redshift 
SNe that are needed to construct a cosmologically inter- 
esting Hubble diagram. Given the unlikely prospects for 
significantly reducing the Hubble scatter, we take a dif- 
ferent approach here and explore models to describe the 
scatter in more detail. Such models of intrinsic scatter 
can be used to evaluate and constrain systematic uncer- 
tainties from assuming an incorrect model, and possibly 

11 Within the SALT-II model framework, these three parameters 
are mj , x\ and c. 



lead to a better understanding of the underlying wave- 
length dependence of SN brightness variations. 

In this work we demonstrate a method for evaluating 
models of intrinsic scatter by computing three scatter- 
dependent dispersion variables and making the follow- 
ing data/MC comparisons: 1) the traditional Hubblc- 
diagram residual, 2) B — V — c, where B — V is the true 
rest-frame color and c is the fitted color parameter from 
the SALT-II model, and 3) the photo- z residual. In the 
hypothetical limit of observing SNe la with infinite pho- 
ton statistics and no intrinsic scatter, the distribution 
for each variable would be a Dirac-(S function. Simula- 
tions that include fluctuations from photon statistics, but 
no intrinsic-scatter model, underestimate the measured 
dispersion in these variables. A viable model of intrin- 
sic scatter must predict the dispersion for each variable 
and for multiple data sets. These three variables do not 
constitute an exhaustive list of photometric observables; 
for example one could examine other rest-frame colors 
(U — B, V — R), correlations among colors, and the de- 
pendence of the scatter on redshift, stretch, and color. 
With our current statistics and signal-to-noise we limit 
this initial study to the three variables described above, 
but larger and higher-quality SN samples from current 
and future surveys should enable a more thorough study. 

Three classes of wavelength-dependent intrinsic scatter 
models are investigated. First we try purely phenomcno- 
logical functions of rest-frame wavelength with param- 
eters tuned to match observations. The second class 
is based on measurements from data. The third class 
uses theoretical explosion models (KRW09) to perturb 
the SALT-n spectral model. 

This work is part of the SDSS+SNLS joint analysis, 
and the data sets used here include 251 spectroscopically 
confirmed SNe la from the t hree-year Sloan Digit al Sky 
Survey-II sample (SDSS-II: iFrieman et all 120081 ). and 
another 191 spectroscopically confir med SNe la from th e 
Supernova Legacy Survey fSNLS3: iConlev et al.ll20Tll) . 
All simulations and lig ht curve fitting are done with 
SNANA version vl0_07 (jKessler et al.1 l2009bl ) and the 
SALT-ii light curve model (G10). 

The outline of the paper is as follows. The data sam- 
ples are described in and the simulation and intrinsic 
scatter models are described in The determination 
of each dispersion variable is in SJH and the resulting 
data/MC comparisons are in Sj5] along with some sys- 
tematics tests. Finally, in ^Hlwe investigate the potential 
Hubble diagram bias from using an incorrect model of 
intrinsic scatter. 

2. THE SDSS-II AND SNLS DATA SAMPLES 

We use two SN la data samples that are well cali- 
brated with ~ 1% photometric precision, and that span 
complementary redshift ranges. The lower redshift SNc 
(z < 0.4) are from th e full 3-season SDSS-II sam- 
ple (|Frieman et al.l 120081) . and the higher redshift SNe 
(z < 1) are from the the pub licly available 3-year SNLS3 
sample (jConlev et al.ll201ll) . Below we give a brief de- 
scription of these samples. 

Th e SDSS-II Supern ova Survey used the SDSS cam- 
era tiGunn et all 11998ft on the SPSS 2.5-meter tele- 
scope ()Gunn et al ] |200a lYork et alj [2000]) at the Apache 
Point Observatory to search for SNe in the Fall seasons 
(1 September through 30 November) of 2005-2007. This 
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survey scanned a region (designated stripe 82) centered 
on the celestial equator in the Southern Galactic hemi- 
sphere that is 2.5° wide and runs between right ascen- 
sions of 20 hr and 4 hr , covering a total area of 300 sq. deg 
with a typical cadence of every four nights per region. 
Images were obtaine d in five broad passbands, ugriz 
(jFukugita et alJl996T) . with 55 s exposures and processed 
throu gh the PHOTO photometric pipeline (|Lupton et al.l 
120011) . Within 24 hours of collecting the data, the images 
were searched for SN candidates that were selected for 
spectroscopic observations in a program involving about 
a dozen telescopes. The SDSS-II Supernova Survey dis- 
covered and spectroscopically confirmed a total of ~ 500 
ty pe la SNe. Det a ils of the S DSS-II SN Survey are given 
in lFrieman et al.l (|2008ft and iSako et~"aH (|2008l ). and the 
procedures for spectroscopic id entification and red shift 
determinations arc described in lZheng etall (|2008t ). 

The SN photometry for SDSS-II is based on Scene 
Mode l Photometry (SMP) described in iHoltzman et al.l 
(|2008j ). The basic approach of SMP is to simultaneously 
model the ensemble of survey images covering a SN loca- 
tion as a time-varying point source (the SN) and sky 
background plus time- independent galaxy background 
and nearby calibration stars, all convolved with a time- 
varying PSF. The fitted parameters are supernova posi- 
tion, supernova flux for each epoch and passband, and 
the host galaxy intensity distribution in each passband. 
The galaxy model for each passband is a 20 x 20 grid 
(with a grid-scale set by the CCD pixel scale, 0.4" x 0.4") 
in sky coordinates, and each of the 400 x 5 = 2000 galaxy 
intensities is an independent fit parameter. As there is 
no pixel re-sampling or image convolution, the procedure 
yields correct statistical error estimates. 

The SNLS was a 5-year survey covering four 1 sq. deg 
fields using the MegaCam imager on the 3.6-meter 
Canada- France-Hawaii Telescope (CFHT). Images were 
taken in four bands similar to those used by the SDSS: 
9m j^MiU/i z m j where the subscript M denotes the 
MegaCam system. The SNLS exposures were ~ 1 hr 
in order to discover SNe at redshifts up to z ~ 1. The 
SNLS images were processed in a fashion similar to the 
SDSS-II so that spectroscopic observations could be used 
to confirm the identities and determine the redshifts of 
the SN candidates. Additional informa tion about the 
SNLS can be found in lAstier et al.l ([2006D and references 
within. 

The SNLS3 SN photometry is based on a simultane- 
ous fit of the SN flux and position, a residual sky back- 
ground per image, and a galaxy intensity map. Images 
are resampled to the same reference pixel grid prior to 
the fit. The SN+galaxy image model is PSF- matched to 
the resampled images. Only sky noise is included in the 
photometric uncertainties (host galaxy and source noise 
are negligible for most supernovac). Because resampling 
introduces pixel correlations, the uncertainties ignoring 
correlations are scaled such that the reduced x 2 is one 
when assuming a constant SN flux per night. 

To ensure good quality fits to the light curves, the fol- 
lowing selection criteria are applied to both the SDSS-II 
and SNLS3 data samples, 

• At least one epoch before the epoch of peak bright- 
ness in the B band (defined as T res t = 0) 

• At least one epoch with T rcst > 10 days. 



• At least three filters with an observation that has 
a signal-to-noise ratio (S/N) above 8 

• At least 5 observations in the fitted epoch range 
— 12 < T r est < +25 days. The maximum T rost is 
set by the range for one of the models (KRW09) of 
intrinsic brightness variation. 

• After fitting each light curve to the SALT-n model 
(' 14.11) . we require the SNIa fit probability to be 
Pfit > 0.02, where Pfn is computed from the fit-% 2 
and the number of degrees of freedom. 

• Color excess from M ilky Way Galactic extinction 
(|Schlegel et al]ll998h is E(B - V) < 0.2. 

The sample statistics after these requirements are given 
in g~3 

3. SIMULATIONS 

We use the SNANA MC code (jKessler et al.l l2009bf ) to 
generate realistic SN la light curves that are analyzed in 
exactly the same manner as the data. The MC is used 
to make detailed comparisons with the data using differ- 
ent models of intrinsic SN la brightness variations. All 
simulations are based on a standard ACDM cosmology 
with w = — 1, Om = 0. 3, f^A = 0-7. Details of the sim- 
ul ation are desc r ibed inlKessler et ah (2009b) and in §6 
of iKessler etaLl (|2009al hereafter K09); here we give a 
brief overview. 

Simulations are generated using the SALT-n model 
(G10) that is based on a time-sequence of rest-frame 
spectra. The spectral model is explained in more detail 
in 2] within the context of light curve fitting. Observer- 
frame magnitudes are computed by redshifting the rest- 
frame spectrum for each e poch, reddening the spec- 
tra from Galactic extinction (jSchlegel et al.lfl998| ) using 
Rv =3.1, and summing the flux in the appropriate filter- 
response curves. 

To account for non-photometric conditions and varying 
time intervals between observations due to bad weather, 
actual observing conditions are used for both the SDSS- 
II and SNLS surveys. For each simulated observation, 
the noise is determined from the measured point spread 
function (PSF)Q Poisson noise from the source, and sky 
background. Noise from the host-galaxy background is 
included for the SDSS-II simulations where it has a small 
effect at low redshifts. Additional details of the simula- 
tion of noise are given in §3.11 The simulated flux in 
CCD counts is based on a mag-to-flux zcropoint, and a 
random fluctuation drawn from the noise estimate. 

The parent distributions of the SALT-II stretch (x\) and 
color (c) are well described by an asymmetric Gaussian 
that is a function of three parameters, 

e [-( Xl - Sl f/2a*_] Xx < - t (1) 

e [-(*i-xO a /a^l Xx > xt (2) 

and a similar function with x\ — > c. The parameters for 
each distribution are shown in Table Q] After account- 
ing for Malmquist bias, we find that the higher-redshift 
SNLS3 sample is slightly brighter and bluer compared 

12 The PSF is described by a double-Gaussian function. 
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to the SDSS-II sample. This difference is expected from 
previous results showing that younger star-former galax- 
ies host brighter /bluer SNe la than older passive galax- 
ies (jSullivan et al.ll2006tlSmith et alJl2012f >. The younger 
star-forming galaxies are more abundant at higher red- 
shifts, thus qualitatively explaining the brightness differ- 
ence between the two surveys. While the variation in the 
color population may be somewhat speculative, we show 
in the systematics analysis f i)5.1[) that the simulated in- 
trinsic scatter is rather insensitive to the parameters de- 
scribing the parent populations in Table [J 



TABLE 1 

Asymmetric Gaussian parameters to describe the parent 
distribution of x\ and c. 





value at 






Generation 




peak prob 




°+ 


Range 


xi (SDSS-II) 


0.5 


1.4 


0.7 


-5,-1-3 


xi (SNLS3) 


0.5 


1.3 


0.7 


-5,-1-3 


c (SDSS-II) 


0.0 


0.08 


0.13 


-0.4, +0.6 


c (SNLS3) 


-0.04 


0.06 


0.14 


-0.4, +0.6 



The simulation includes a detailed treatment of the 
search efficiency, including spectroscopic selection effects. 
For the SDSS-II, the search-pipeline efficiency has been 
measured separately for each q,r,i filter usi ng fake SNe 
inserted into the images (jDildav et all 120081) . The spec- 
troscopic selection efficiency (e sp0 c) has been estimated 
from matching data/MC distributions for redshift and 
for the fitted observer-frame magnitudes at the epoch of 
peak brightness. e S pcc is adequately described as a func- 
tion of peak r-band magnitude and the peak color g — r. 
These efficiency functions are available in tabular formF^I 
For the SNLS3, e spcc has been evaluated as a function of 
eak i m band magnitude (Mj) in Fig. 9 of iPerrett et al.l 
2010D . For the SNANA simulation we parameterize this 
function as 



e spoc S3 — ^ 0-5 



1 



■ tan 



24.3 - Mi 
02 



xC e (Mi) (3) 



where C e (Mj) = 1 for Mj < 23 and C e (Mj) = exp[(23 - 
Mj)/0.63] for Mi > 23. The function in parentheses is a 
first-order estimate and (7 e (Mj) is a correction obtained 
from a fit to the data/MC ratio as a function of Mj. 

For this analysis we generate MC samples with sizes 
corresponding to six times the data statistics. The qual- 
ity of the simulation for each sample is illustrated with 
several data/MC comparisons in Figs Q][2J the overall 
agreement is good. 

3.1. Simulation of Noise 

Since the three scatter-dependent variables (Hubble 
scatter, B — V — c photo- z precision) are sensitive to the 
flux uncertainties, it is important to accurately simulate 
these uncertainties. The simulation strategy is to first 
calculate the uncertainties from a model based on mea- 
surements of the sky level and PSF. To accurately check 
the model, the true uncertainty for each observation in 
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Fig. 1. — Comparison of distributions for SDSS-II data (dots) 
and MC (histogram), where each MC distribution is scaled to have 
the same sample size as the data. The distributions are redshift, 
number of degrees of freedom in the SALT-II light curve fit, maxi- 
mum rest-frame time difference (gap) between observations, maxi- 
mum S/N, fitted SALT-II color (c) and stretch parameter (xi). The 
bottom two panels show the mean fitted SALT-II color (c) and shape 
parameter (a;i) versus redshift. 
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Fig. 2. — Same as Fig[T] except for the SNLS3 sample. 
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the datg(3 is compared to the calculated model uncer- 
tainty. Discrepancies between the true and calculated 
uncertainties are corrected by fitting for ad-hoc parame- 
ters. The simulated uncertainty model in photoelectrons 
(csim) is given by 

°SiM = F + (A ■ b) + (qF) 2 + ((To • io°' 4 ' ZPTpc ) 2 

+ °host ( 4 ) 

where F is the flux, A = [2tt / PSF 2 (r, eydr]- 1 is the 
noise-equivalent area, b is the effective sky level including 
dark current and readout noise, and q and oo are fitted 
ad-hoc parameters. ZPT pe is defined such that the num- 
ber of CCD photoelectrons for a point source of magni- 
tude m is given by io~°- 4 ( m ~ ZPT P'=); thus the <jq term is 
independent of the PSF, sky level and host-galaxy. The 
quantity cthost is simulated for the SDSS-II sample us- 
ing a library of galaxies that have a spectroscopic redshift 
and a well measured profile consistent with either expo- 
nential or de Vaucouleurs. 

To check the uncertainty calculation, ctsim is computed 
for each epoch in the data and compared to the measured 
uncertainty ctdata- The left panels in Fig. [3] show that 
the first two terms, F + Ab, are not adequate to repro- 
duce the observations. The right panels in Fig. [3] show 
that the fitted oo term results in good agreement over a 
wide range of PSF values. A separate oo value is evalu- 
ated for each filter and for each sample (Table [5]) . The 
quadratic term q is sensitive to large flux values with 
S/N ~ 10 2 . The value of q is obtained from minimiz- 
ing x 2 = X^K^Sim/o-data);; - I] 2 , where the sum (s) is 
over log 10 (S/N) bins; q ~ 0.01 for the SDSS bands, and 
q ~ 0.001 for the SNLS bands. For SNLS, it is difficult to 
interpret this low value of q because uncertainties on the 
SN flux (Poisson noise and flat-fielding noise) only arise 
via the normalization of errors based on the intra-night 
flux scatter. 

Finally, note that the terms F + (A ■ b) + Ohqst are 
determined from observations and first principles, while 
q and oo are empirically determined parameters. The 
q term corresponds to a zeropoint uncertainty. The oo 
term is not understood, although this term works sur- 
prisingly well for both the SDSS-II and SNLS3 surveys 
even though the respective photometry codes are inde- 
pendent. 



TABLE 2 

Fitted <tq values for SDSS-II and SNLS3. 



Filter 


10 11 x 
SDSS-II 


do for 
SNLS3 


u 


28 




9 


10 


0.22 


r 


15 


0.30 


i 


23 


1.08 


z 


63 


1.62 



True data uncertainties are from the photometric pipeline 
codes describe in § [2] 
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Fig. 3. — Left panels show the r-band flux-uncertainty distribu- 
tion for the data (black dots) and for the calculation (histogram), 
CT SIM = VF + Ab with (To — and q = 0. Each descending plot 
is for a different PSF-range (FHWM, arcsec) as indicated on the 
plot. Right panels show the ctsim calculation including the best-fit 
(To value labeled on the plot and shown by the arrow. The flux 
unit shown is that used for light curve fits (not photoelectrons). 

3.2. Intrinsic Scatter Models 

The intrinsic scatter models are summarized in Tablc|3l 
These models are defined as wavelength-dependent per- 
turbations to the SALT-n spectral model, and these per- 
turbations average to zero so that the underlying SALT- 
II model is not changed. All models are independent 
of redshift, and only the explosion models from KRW09 
depend on epoch. We begin with the phenomenological 
functions (see "FUN" prefix) with parameters arbitrar- 
ily chosen to increase the scatter. The coherent model 
(FUN-COH) assigns a single magnitude shift for all wave- 
lengths; for each SN this shift is given by a Gaussian- 
random number with ctcoh = 0.13 mag. The other 
two FUN functions are designed to probe a wider va- 
riety of wavelength-dependent scatter with a coherence 
length of a few hundred A. First a sequence of nodes is 
defined at 1000 A intervals in the rest frame. An in- 
dependent Gaussian random scatter is selected at each 
node with er no dc = 05500 cxp[— (A no dc ~ 5500)/3000] so 
that there is more scatter at bluer wavelengths. The 
variation is the same at all epochs. A continuous func- 
tion of wavelength is constructed by connecting the node 
values with sine functions so that the derivative is zero 
at each node. FUN-COLOR is defined with 05500 = 0.06 
and is shown in Fig.g]for a few simulated SN. FUN-MIX 
is defined with 05500 = 0.045 along with a coherent term 
o-coh = 0.09 mag. 

The next two models (G10 and Cll) are based on mea- 
surements from data combined with assumptions needed 
to create a model that is a continuous function of wave- 
length. The G10 error model was obtained as part of the 
SALT-n training process in which they minimized the like- 
lihood of light curve amplitude residuals using a paramet- 
ric function of central rest-frame wavelength, assuming 
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TABLE 3 

Summary of intrinsic scatter models. A ros t refers to 

REST-FRAME WAVELENGTH. 



iviociei 






description 


NONE 


Poisson noise only 


FUN-COH 


coherent variation at all A r0 st 




(°-COH = 0.13 mag) 


FUN-COLOR 


^nodo = 0.06exp[-(A node - 5500)/3000] 




(Fig. m 


FUN-MIX 


<r nodo = 0.045 exp[-(A nodc - 5500)/3000] 




and (jcoh = 0.09 


G10 


O-COH = °' 09 ma § and °node 




from solid curve in Fig. 8 of G10 


C11.0 


UBVRI correlation matrix, pui y = 0. 


C11.1 


same for UBVRI, pjji jj = +1. 


C11.2 


same for UBVRI, piji [/ = —!. 


(KRW09 2D Explosion Modcls a ) 


iso3_dcl 


Mgnit = 60, Ka min = 1 


iso6_dcl 


Mgnit = 100, Ka min = 1 


iso6_dc2 


Mgnit = 100, Ka min = 250 


iso6_dc3 


Mgnit = 100, Ka min = 750 


iso6_dc4 


AT ignit = 100, Ka min = 1500 


iso6_dc5 


jV ignit = 100, Ka min = 2250 


iso8_dc3 


^ignit = 150, Ka min = 750 


asyml_dc3 b 


^ignit = 120, Ka min = 750 



a "iso" = isotropic and "dc" = detonation criteria. 

b This mis-labeled asym model is really an isotropic model. 



uncorrelatcd residuals in different passbands. The result- 
ing wavelength-dependent function (Fig. 8 of G10) has 
approximate values of 0.07, 0.03, 0.02, 0.03, 0.06 mag at 
the U, B, V, R, I central wavelengths, respectively. This 
function is not intended to represent a wavelength de- 
pendent scatter, but rather it is a model of independent 
broadband scatter as a function of central wavelength. 
To translate this broadband model into a wavelength 
model, independent random scatter values (cr no de) are se- 
lected every 800 A, and these node values are connected 
with the same sine-interpolation that is used for the phe- 
nomenological "FUN" functions. Since this procedure re- 
duces the resulting broadband scatter, er no de i s multiplied 
by 1 + (A rest - 2157) /9259 so that the simulated UBVRI 
broadband scatter matches the G10 function. In addi- 
tion to a wavelength-dependent function, the G10 model 
includes a coher ent term, aeon = 0.09 mag. 

The model of IChotardl (l201ll) : IChotard et all (l201lL 
hereafter Cll) is based on a covariance scatter ma- 
trix among the UBVRI filter passbands, and is derived 

i_ 0.3 q 1 t-, 1 t 1 q 1 




2500 5000 7500 2500 5000 7500 2500 5000 7500 2500 5000 7500 

X (A) X (A) X (A) X (A) 



Fig. 4.— Illustration of FUN-COLOR for four simulated SNe. 
Black lines show the A res t-dependent variation. Solid dots show the 
Gaussian-random variations at the nodes, and dashed lines show 
±1 sigma vs. A rcst . 



from an analysis of spectral correlations using high qual- 
ity sp ectra from the Supernova Factory (|Aldering et al.l 
120021) . The broadband covariance model is translated 
into a wavelength-dependent model as follows. First, 
the model is extrapolated to wavelengths below the U 
band (3600 A) by defining an ad-hoc U' filter with cen- 
tral wavelength A b s = 2500 A. The G10 scatter value of 
Cnodo = 0.59 mag is used for U', and we model three dif- 
ferent assumptions for the reduced correlation between 
U' and U: pw.u = (incoherence 11_0), pw.u = +1 
(Cll_l), and pu'.u = — 1 (Cll_2). For each simulated 
SN, six random magnitude shifts are selected according 
to the Cll correlation matrix in upper half of Table SJ 
these shifts are assigned to the central U'UBVRI wave- 
lengths. A continuous function of wavelength is obtained 
by interpolating these six points with a sine function, 
similar to the FUN-COLOR interpolation in Fig. g] Fi- 
nally, the scatter function is multiplied by 1.3 to com- 
pensate for the fact that the wavelength interpolation 
reduces the broadband covariances. The correlation ma- 
trix realized by the simulation is shown in bottom half 
of Table|l]for the Cll_l model. The realized correlation 
matrix is slightly different than the input model because 
the input model is described by broadband covariances, 
while the simulated model depends on wavelength. In 
principle a more finely tuned spectral model in the simu- 
lation would result in the exact Cll covariances, but we 
believe that the simple and approximate model used here 
is adequate, especially in light of the large and unknown 
uncertainties on the covariances. 



TABLE 4 

Reduced correlation matrix from the Cll_l model 51 and 

REALIZED FROM THE SIMULATION. 





U' 


U 


B 


V 


R 


I 








(From 


Cll) 






u' 


+1.00 


+ 1.00 


-0.12 


-0.77 


-0.91 


-0.22 


u 


+1.00 


+ 1.00 


-0.12 


-0.77 


-0.91 


-0.22 


B 


-0.12 


-0.12 


+ 1.00 


+0.57 


-0.24 


-0.89 


V 


-0.77 


-0.77 


+0.57 


+1.00 


+0.53 


-0.40 


R 


-0.91 


-0.91 


-0.24 


+0.53 


+1.00 


+0.49 


I 


-0.22 


-0.22 


-0.89 


-0.40 


+0.49 


+1.00 


Diag b 


0.59 


0.06 


0.04 


0.05 


0.04 


0.08 






(realized in 


simulation) 




V 


+1.00 


+0.99 


+0.12 


-0.78 


-0.96 


-0.28 


u 


+0.99 


+ 1.00 


+0.20 


-0.73 


-0.96 


-0.34 


B 


+0.12 


+0.20 


+ 1.00 


+0.40 


-0.27 


-0.93 


V 


-0.78 


-0.73 


+0.40 


+1.00 


+0.72 


-0.34 


R 


-0.96 


-0.96 


-0.27 


+0.72 


+1.00 


+0.36 


I 


-0.28 


-0.34 


-0.93 


-0.34 


+0.36 


+1.00 


Diag a 


0.59 


0.073 


0.036 


0.051 


0.042 


0.083 



a See Table 14.2 oflChotardl pjll PhD thesis). 
b Diag = VCOVii with i = U', U, B, V, R, I are from the PhD 
thesis and differ slightly from those given in Chotard et al.l l|2011I V 

The final class of brightness variations is based on 2D 
explosion models with random ignition points (KRW09), 
followed by radia tive transfer calcu lations using the 
SED0NA program (lKasen et al.l 120061 ). Isotropic mod- 
els are obtained from ignition points that are randomly 
placed throughout the white dwarf (WD), while asym- 
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metric models are obtained from ignition points within 
a cone whose apex is at the center of the WD. Both 
the isotropic and asymmetric models result in explosion 
asymmetries and a viewing angle dependence that con- 
tributes significantly to the intrinsic scatter. The width- 
luminosity relation is related to the number of ignition 
points (-/Vjgnit) because iVj gn it affects the amount of pre- 
expansion before detonation, and hence the amount of 
56 Ni produced in th e explosion. In a recent study by 
iBlondin et all (j201lL hereafter Bll), detailed compar- 
isons between data and the KRW09 models were made. 
They conclude that the KRW09 models with the best 
spectroscopic agreement also have the best photometric 
agreement, and they identified a subset of eight models 
with the best agreement to data. Here we use these same 
eight models; they are shown in Table |3] along with a 
few parameters describing the number of ignition points 
and the detonation criteria. All of these models have 
an isotropic distribution of ignition points, and Bll note 
that radial fluctuations in isotropic models can lead to 
significant viewing angle asymmetries. 

We initially used these KRW09 models in the SNANA 
simulation to generate light curves corresponding to the 
SDSS-II and SNLS3. While the data/MC comparisons 
are visually impressive, the simulated light curves are not 
adequate for this study because the SALT-n light curve 
fits are in general rather poor. This trend of poor light 
curve fits was also noted in Bll. 

Instead of attempting an absolute prediction with the 
KRW09 models, we have instead used these models as a 
perturbation on the SALT-n model. In short, the SALT-n 
model describes the stretch and color relations, while the 
KRW09 models are used to describe the intrinsic scatter. 
The spectral flux (F) is given by 



F — -FgALT2 x 



Ft 



KRW09 



(random cost 



^KRW09 (COS <f> = 0) 



(5) 



where <fi is the viewing angle. The corresponding mag- 
shifts are illustrated in Fig.[5]as a function of wavelength 
for a two extreme viewing angles. 

4. ANALYSIS 

Here we describe the determination of the three 
scatter-dependent quantities used to evaluate models of 
intrinsic brightness variations. All analyses are based on 
light curve fits using the SALT-n model. 

4.1. Review o/salt-ii Model 

The SALT-n SN la model flux is a function of wave- 
length (A) and time (<) in the rest-frame, 

F(t, A) = .t [M (t, A) + xtMxit, A)] x exp[c-OL(A)] (6) 

where the spectral sequences (Mq and Mi) and color 
law (CL(A)) are derived from the training in G10. Syn- 
thetic photometry in the observer frame is obtained by 
redshifting Eq. [6] and multiplying by the filter response 
and Galactic transmission. The overall scale (xo), stretch 
(xi), color (c) and time of peak brightness (to) are deter- 
mined for each SN in a light curve fit that minimizes a y 2 
based on the difference between the data and synthetic 
photometry. Eq. M is valid for 2000 < A < 9200 A, and 
the model is valid for observer-frame filters that satisfy 



n3 
E 

O 
EC 



-0.5 

0.5 


-0.5 



z iso3 dd 


iso6 dd 


i iso6 dc2 


iso6 dc3 


\ iso6 dc4 


iso6 dc5 


i i i i i i i i i i i 


V asyml dc3 



2000 



6000 80002000 

MA) 



6000 8000 

MA) 



Fig. 5. — Magnitude shift applied to the SALT-II spectral model 
vs. wavelength for the KRW09 models. The solid and dashed 
curves reflect different viewing angles as indicated in the legend 
above the plots. 



where A Q b s is the central wavelength of the filter. 

An effective -B-band magnitude is defined to be = 
— 2.5 log 10 (xo) + 10.635; this is the observed magnitude 
through an idealized filter that corresponds to the B 
band in the rest-frame of the SN. The fitted distance 
modulus is given by 



//fit = rriB — M + axi — (3c 



(8) 



where a, f3 and M are determined from a global fit to 
all of the SNe using the "SALT2mu" program described in 
Mil and below in $6] 

4.2. Hubble Scatter 

The well known Hubble scatter is defined as the disper- 
sion on A/z, the difference between the fitted (measured) 
distance modulus and the distance modulus calculated 
from the best-fit cosmological parameters. To simplify 
the analysis here we don't fit for the a and (3 parame- 
ters, nor do we fit for the best-fit cosmological parame- 
ters. Instead we compute the dispersion of 



A/i = yUfi t (Eql5]) - ^calc(z, ^M, ^A, w) 



(9) 



2800 < A obs /(l +z) < 7000 



(7) 



where a = 0.11, /3 = 3.2, M = -19.36 (H = 
70 km/sec/Mpc), and /x ca ic is the calculated distance 
modulus assuming a ACDM cosmology with w = — 1, 
SIm = 0.3, £1a = 0.7. Although the fitted a and (3 may 
have given slightly different AyU values, the resulting bias 
is more than an order of magnitude smaller than the dis- 
persion, and hence the impact of this approximation is 
negligible. 

4.3. Color Precision 

The color precision test compares the fitted SALT-ii 
color (c) to the true B — V rest-frame color at the epoch 
of peak brightness. Since the fitted color is really a color 
excess, c = E(B — V), and the color also depends slightly 
on the stretch, the true B — V color does not exactly 
correspond to c. A numerical examination of the model 
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shows that with no intrinsic scatter, 

B — V = c' = (1.016 x c) + (xi/1250) + 0.0232 (10) 

and therefore we examine the dispersion on B — V — c' . 
The dispersion on c — c' is ~ 0.001, more than an order of 
magnitude smaller than the dispersion on B — V— c' . The 
evaluation of c' is from simply plugging the fitted c and 
x\ values into Eq. 1101 The naive rest-frame magnitudes 
Mg and M v are obtained from Eq. [6] using the best-fit 
parameters (c, xq, x\) and the B and V filter transmis- 
sion functions. However, these naive magnitudes are not 
necessarily the true values if there are intrinsic color vari- 
ations. To obtain a better approximation for the magni- 
tudes we fit only the two nearest observer-frame bands 
that bracket the B or V band in wavelength. 

The details of the fitting procedure are as follows. First 
a normal fit is done using all filters to determine the fit 
parameters (to-,c,Xo, x\). For each rest- frame band one 
additional fit is performed using only the two nearest 
observer- frame bands and holding to and x\ fixed from 
the normal fit. The floated color (c) and distance (xq) 
parameters provide the flexibility to fit both observer- 
frame bands regardless of how much intrinsic color vari- 
ation exists. The two-band fit parameters are c B , x B for 
the B band, and c v ,x^ for the V band. After finishing 
both two-band fits the B — V color is computed as 

B-V = M*(T B , t , xi = 0, c B , x B ) 

-M*(T V , t , Xl =0, c v , x%) (11) 

where M* is the magnitude computed from Eq. [6] using 
filter-transmission functions Tb,v, and with x\ = so 
that all B — V colors correspond to a SN la with the 
same stretch. 

This fitting procedure was tested on an SNLS3 sim- 
ulation in which the maximum S/N was artificially set 
to 1000 for every SN regardless of redshift. The rms on 
B — V — c' is 0.002 mag, an order of magnitude smaller 
than the observed dispersion. 

4.4. Photoz Precision 

The photo-z precision is based on the difference be- 
tween the SN redshift determined from broad-band pho- 
tometry and the more precise spectroscopic redshift. The 
basic photo-z method is to extend the usual methods of 
fitting light curves to include the redshift as a 5 th fit 
parameter. Particular attention is needed to estimate 
initial parameter values near those corresponding to the 
global minimum-;^ 2 , and to iteratively determine which 
filters satisfy Eq. [7] Details of the photo-z fitting process 
are given in K10. 

There are two modifications in our photo-z fitting pro- 
cedure compared to K10. The first change is that we use 
the known spectroscopic redshift as the initial estimate 
in order to reduce catastrophic outliers. The fitting task 
has thus been changed to find a local photo-z minimum 
near the true redshift instead of searching the entire red- 
shift range for a global minimum. The second change is 
related to estimating the initial parameter xq for each 
color value along the coarse-grid search in color. In K10, 
xo at each grid point was calculated using the current 
color, photo-z and a reference cosmology. Here we ana- 
lytically minimize for xq, making the fitted photo-z less 
sensitive to the absolute brightness. To check that the 



fitted photo-z depends only on the SN colors we have ap- 
plied this method to simulations with no intrinsic scatter 
and with the coherent scatter model (see COH entry in 
Table EJ); the photo -z precision is the same in both cases. 

4.5. Statistics Summary 

After applying the selection requirements in Sj2j along 
with the light curve fitting requirements for each disper- 
sion variable, the number of SNe la for each sample and 
for each dispersion variable is shown in Table [5] 

The smaller photo-z samples arise from a light curve 
fitting requirement. For each successive fit iteration, 
observer-frame filters are added or dropped based on 
which filters satisfy the SALT-n wavelength range in Eq.[7] 
with z = photo-z. If any filter fails this wavelength re- 
quirement after the last fit iteration, the SN is rejected; 
this requirement avoids fitting to wavelength regions in 
which the SALT-ii model may be poorly defined. 

The smaller SNLS3 sample for the B — V — c' analysis is 
due to SNe la at z > 0.7; for these objects the observer- 
frame i and z bands no longer bound the rest-frame V 
band. 



TABLE 5 

Number of SNe Ia After Selection Requirements. 



Survey 


Number of SNe Ia for: 
Hubble resid B — V — c' photo-z 


SDSS-II 
SNLS3 


251 245 208 
191 120 176 


4.6. 


Quantifying the Dispersions 



The data and MC dispersions are measured from the 
following variables, 

A M = (/ifit - /icalc)/ (1 + Z) 

A c = (B-V-c')/(l + z) 

A z = {Z phot - Z spcc )/(l + z) (12) 

where the (1 + z) _1 factor is included to reduce the 
redshift-dependent variation from measurement uncer- 
tainties. These three quantities are shown in Fig. [6] for 
the data and MC, and for both the SDSS-II and SNLS3 
samples. The MC includes only Poisson noise (no intrin- 
sic variation), and hence the data-MC difference in the 
width illustrates the size of the intrinsic component that 
is needed. The A M comparison shows the most obvious 
discrepancy. The A z and A c discrepancies are more sub- 
tle, indicating that the effect of color variations is smaller 
than the coherent variation. 

To quantify the dispersions we compute the median, 
Ma = median | Aa, |, where x = /i, z,c indicates the vari- 
able type. In particular, we compute the MC/data ratio 
of medians, 

i?MC/Data = M^/M^ . (13) 

With the correct model of brightness variations we ex- 
pect i?MC/Data = 1 f° r au three variables and for both 
surveys. 
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Fig. 6. — Comparison of A M , A c and A z distributions for data 
(dots) and MC (histogram). The upper plots are for the SDSS-II; 
lower plots are for SNLS3. The MC includes Poisson noise but no 
intrinsic brightness variation, and each MC distribution is scaled 
to have the same sample size as the data. 



The uncertainty on the median is calculated as follows. 
For N supernovae, the statistical uncertainty on N/2 is 
>~N /2. The median uncertainty (ctm) is defined 



O {N/2} 

such that N — 



C(7V/2) 



values of | Ax | lie below M — a M and 



N +<J(n/2) values lie below M + a^- For a rapidly falling 
distribution we typically find that cr^ > er^. Here we 



define a symmetric uncertainty, cm = (erj^ + cr^ I )/2. 

As a numerical crosscheck we analyze SDSS-II simula- 
tions in which the exposure time is adjusted for each SN 
so that the maximum S/N is 10 4 . The resulting disper- 
sions, defined simply as a Gaussian fitted a, are 0.0003, 
0.002, and 0.001 mag for the three variables (Eq. [T2"J) . 
respectively; these dispersions are more than an order of 
magnitude smaller than the dispersions observed in the 
data. 

5. RESULTS 

The MC/data ratio of medians, i?MC/Data (Eq. fT5|) . is 
shown in Fig. [7] for all of the models in Table [31 and 
for both surveys. With no model of intrinsic scatter, 
i?MC/Data is well below unity in all cases. Adding a coher- 
ent scatter (FUN-COH) predicts the Hubble dispersion 
(A^), but has no impact on the color and photo-z disper- 
sion. The FUN-COLOR model almost predicts the Hub- 
ble dispersion, but may overestimate the color precision. 
FUN-MIX has been artificially tuned to predict the dis- 
persion in all quantities, although the photo-z dispersion 
may still be underestimated. The G10, C11_0 and Cll_l 
models provide decent predictions, with a slight underes- 
timate in the photo-z dispersion. The Cll_2 model un- 
derestimates the Hubble dispersion. Recall that the G10 
model includes only positive correlations, mainly from 
the coherent term ctcoh = 0.09, while the Cll model in- 
cludes both positive and negative correlations. This G10 



versus Cll comparison illustrates that there can be sig- 
nificant degeneracies among models of intrinsic bright- 
ness variations. The KRW09 models give a poorer de- 
scription of the dispersion because the Hubble dispersion 
is always underestimated. 

5.1. Systematics Tests 

Here we describe some systematics tests to demon- 
strate the robustness of the results in Fig. [7J We use 
the FUN-COH scatter model as the reference simulation 
for these tests which are summarized in Fig. [5] For each 
test a change is applied to the simulation and then ana- 
lyzed in exactly the same manner. 

The first test is based on the precision in the flux un- 
certainties in the data. For pre-explosion epochs in which 
the true SN flux is known to be zero, examining the S/N 
distribution shows that the uncertainties are accurate to 
within 5%. The test labeled ctsim x 1-05 corresponds to 
a simulation with 5% larger uncertainties on all of the 
fluxes. 

The next set of tests is based on a 0.02 mag zeropoint 
change in each filter (SZPT gr i Z ). Note that this change 
is two times larger than the uncertainty reported by each 
survey team. 

Uncertainties on the Galactic extinction are examined 
by first increasing the estimated 16% scatter to 24% 
(xl.5 MW-Gal), and then increasing the reddening pa- 
rameter by 10%, to Ry — 3.4. 

The next test is based on changing j3 from 3.2 to 2.5, 
a 5(7 change from G10. 

The next two tests are based on changing the popula- 
tion parameters for x\ and c (Tablc[T]). The simulated x\ 
population is shifted toward faster-declining light curves 
by setting cr + = 0.5 and <7_ = 1.8 (compare to nomi- 
nal parameters in Table [T]). The simulated color pop- 
ulation is shifted toward the red by setting er + = 0.18 
and cr_ = 0.05. For these systematic tests, the resulting 
data/MC comparisons for stretch and color are shown in 
Fig. El the data and MC are clearly discrepant. 

For all of these systematic tests, -RMC/Data remains sig- 
nificantly below unity for the B — V — c' and photo-z 
variables. We also note that the independent SDSS-II 
and SNLS3 results are consistent, showing consistency 
over different redshift ranges. 

6. IMPACT OF INTRINSIC SCATTER MODEL ON 
THE HUBBLE DIAGRAM 

Here we investigate the potential Hubble diagram bias 
from using an incorrect model of intrinsic scatter. We 
use four intrinsic scatter models that give reasonable 
data/MC agreement in Fig. FUN-MIX, G10, C11.0 
and Cll_l. Recall that data/MC agreement in these dis- 
persion variables means that the intrinsic scatter model 
cannot be ruled out, but the agreement does not ensure 
that the underlying model is correct. 

The Hubble bias is determined from the difference be- 
tween an ideal analysis using the correct intrinsic scatter 
matrix ( £|6.1[) . and a conventional analysis that adds a 
wavelength-independent scatter to bring the reduced x 2 
to unity f t)6.2[) . The ideal analysis is based on simulations 
with the correct model of intrinsic scatter, and thus does 
not reflect a realistic analysis that could be applied to 
data. The conventional analysis, however, reflects a re- 
alistic analysis that has often been applied to data. In 
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for A M ,z, c as indicated in the legend panel 
from Fig. [7J The first panel is for the FUN-COH scatter model 
using the nominal simulation. The systematic test indicated above 
each panel is applied to the simulation using the same FUN-COH 
scatter model. The horizontal dashed lines are the same in each 
panel, and are intended to guide the eye. 



£!6.3l thc Hubble bias is translated into a bias on the dark 
energy equation of state parameter u>, and in £16. 41 the 
biases are reevaluated with Malmquist bias corrections. 

6.1. Determining the Intrinsic Scatter Matrix 

To evaluate the effect of intrinsic scatter in the anal- 
ysis of cosmological parameters, we first need to briefly 
summarize the concept of an intrinsic scatter matrix in- 
troduced in § 2 of Mil. Cosmology fitters in general 
minimize the function 

X 2 = J2 A ^/VstaM 2 + CTinM 2 ) (14) 

i 

where A/i^ is the difference between the fitted and cal- 
culated distance modulus (Eq. [9]) for the i th SN, o- s tat,i 
is the statistical (fitted) error on A^, and <7i nt j is an 
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Fig. 9. — Comparison of fitted stretch (xi) and color (c) distri- 
butions for data (dots) and modified MC (histogram), where each 
MC distribution is scaled to have the same statistics as the data. 
The simulated populations have been modified for systematic tests 
explained in the text. 

ad-hoc parameter defined so that x 2 /-^dof = 1- Since 
fstat.i is computed from a statistical correlation matrix 
between the SALT-ii fit parameters (1715,11, c), Mil pro- 
posed an analogous "intrinsic-scatter covariance matrix" 
(denoted S) to compute o~i n t. Dropping the SN index i, 
the ad-hoc error term is 



'int 



= £0 



-2aSoi-2/3E 0c -2a/3S]i c 



(15) 



w here the subsc ript c orrespondenc e is , 1, c — > rriB,Xi,c, 
= a^ B , = 4T\ and V^Tc = af. All SNIa- 

cosmology analyses to date have used the simplifying as- 



sumption that <7i, 



and ignored the 



other E terms in Eq. [TS] We refer to this method as 
the "ms-only" method, while the E-fit method refers to 
using additional terms in Eq. 1151 The m^-only method 
is valid if the intrinsic scatter is independent of wave- 
length, or if c s tat,i from the light curve fit includes the 
wavelength dependence of the scatter. The FUN-COH 
panel in in Fig. [7] clearly shows that the intrinsic scatter 
cannot be constant (i.e., wavelength independent). Mil 
noted that using the ms-only method can lead to biased 
values of a and f3. Here we go a step further and examine 
biases in simulated Hubble diagrams. 

Since we don't have reliable methods for determining 
E from the data, we determine E from an artificial analy- 
sis using a simulation with the correct model of intrinsic 
scatter but without Poisson fluctuations from the calcu- 
lated measurement uncertainties (Eq. |4]); therefore the 
only source of scatter is from intrinsic variations. Al- 
though Poisson fluctuations are not applied, the uncer- 
tainties are included in the light curve fit-x 2 calculations 
so that the correct filter-dependent weights are used. For 
example, the SDSS u band has relatively poor S/N com- 
pared to the other bands and therefore this passband 
has less weight in determining E. We refer to these sim- 
ulations as "intrinsic-only" to distinguish them from the 
"nominal" simulations that include Poisson fluctuations. 
This intrinsic-only simulation is illustrated in Fig. 1101 for 
the special case with no intrinsic scatter; the simulated 
fluxes lie exactly on the best-fit SALT-n model and they 
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have the correct uncertainties corresponding to real ob- 
servations. 

To better compare the re s ulting bias to the uncertainty 
reported in iSullivan et al.l (|2011| ). the simulations have 
been adjusted to better match the data sample used in 
this SNLS3 analysis. First, the SDSS-II sample size is 
reduced by a factor of three to correspond to the first- 
season sample used in the SNLS3 analysis. The next 
change is that the S/N requirement in the three pass- 
bands (see end of Sj2| is relaxed from 8 to 5. Finally, we 
have included a simulated nearby (z < 0.1) sample as ex- 
plained in the Appendix. To measure biases with good 
precision, the MC sample sizes correspond to 20 times 
the data statistics. 

After performing SALT-n light curve fits on the simu- 
lated intrinsic-only sample, we define = ms(fit) — 
ma (true) and A^ IM = c(fit) — c(true), where "true" in- 
dicates the true value from the simulation and "fit" in- 
dicates the result from a light curve fit. The true values 
are defined by the underlying SALT-ii model before the in- 
trinsic smearing model is applied. The covariance terms 
with the stretch parameter x\ are negligible because the 
intrinsic scatter models are epoch-independent and thus 
do not change the light curve shape; the xi-terms in E 
arc therefore ignored. 





A 



Fig. 10. — Simulated SDSS-II light curve (solid dots) with no 
Poisson noise and no intrinsic scatter. Each curve is the best-fit 
SALT-II model. 

The 2x2 intrinsic scatter matrix is defined to be 

^ B =rms|Aei 
i nt = rms|A? IM | 




E 0c = (A--A 



SIM \ SIM^ 



(16) 

where () indicates the mean value of the enclosed quan- 
tity. Another caveat is that E depends on the redshift 
and on which filters are included in the light curve fit. 
This dependence is linked to the SALT-n color parameter 
(c) that is evaluated by extrapolating a color law to the 
central wavelengths of the B and V passbands. To ad- 
dress this caveat, £ is evaluated as a function of redshift 
and sample as shown in Fig. 1111 A second-order poly- 
nomial function of redshift is adequate to describe the 
components of £. For the FUN-MIX and G10 models, 
a^ B and /3aJ, nt gives a comparable contribution (~ 0.1) 
to cr int . For the Cll models, a"f B ~ 0.03 is much smaller 
than the contribution from /3a, 



hit 



6.2. Fitting for a, j3, and the Distance Moduli 

Using the fitted SALT-n parameters and uncertainties 
from the nominal MC, we use the SALT2mu program 
(Mil) to minimize Eg. 1141 This minimization gives the 
best-fit values of a and /3, along with an independent 
offset (M z ) in each 0.1-wide redshift bin (Eq. |SJ). The 



separate offsets are used to eliminate the /Lt-dependcncc 
on cosmological parameters. These fitted parameters are 
used to compute a distance modulus for each SN, and the 
resulting rcdshift+distance pairs form a Hubble diagram 
that can be fit for cosmological parameters. To properly 
evaluate the fit x 2 using the scatter matrix E, the light 
curve fits have been done with the SALT-n model uncer- 
tainties set to zero (hence reducing <x s t a t,i) because they 
are now included in the SALT2mu fit as part of the o~i n t 
calculation. This modification only affects the SALT2mu 
X 2 and has a negligible impact on the fitted parameters. 

The ms-only fit is based on the traditional method 
of tuning cr™^ such that xV-^dof = 1, while setting all 
of the other covariance terms to zero. The light curve 
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Fig. 11. — Elements of intrinsic scatter matrix (S ) for the 
simulated models indicated above each set of plots. The panels 
show the following parameters vs. redshift: u mB = V^oo (loft), 
3<r c = 3\/^cc (middle), and p mB ,c = ^Oc/(&m B Vc) (right). The 
"int" superscripts have been dropped for clarity. Note that (7m B 
and 3cr c are in roughly the same units as the scatter in the distance 
modulus residual. The simulated SNLS3 (SDSS-II) samples are 
indicated by solid (open) circles. Each solid curves represents a 
second-order polynomial function of redshift used to describe S. 



13 



fits include the SALT-n model errors, and hence the in- 
trinsic uncertainties corresponding to the G10 model are 
included. 

The SALT2mu fit results for both intrinsic scatter meth- 
ods (m^-only and E) are shown in Table [5] For both 
methods the fitted values of a agree well with the sim- 
ulated input, asiM = 0.11. For the FUN-MIX and G10 
scatter models, the fitted values of /3 agree well with the 
simulated input (/?sim = 3.2) for both methods. For 
the Cll models, however, the situation is quite different. 
The ms-only fitted j3 values are too low by 0.3-0.6, and 
the significance of this bias ranges from 6-15 standard 
deviations. The E-fit /3 values are consistent with /?sim- 

The reduced \ 2 values are close to unity for the ttib- 
only method because of the explicit a™^ tuning. How- 
ever, there is no such tuning for the E-fit method and the 
reduced \ 2 are within about 10-20% percent of unity. 
In summary, using the correct intrinsic scatter matrix 
E in the SALT2mu fit results in unbiased /3 values, and 
X /Ndoi — 1- Using the simplistic m^-only method re- 
sults in a significant bias on /3 for the Cll models. 

Since the E-fit method results in unbiased j3 values we 
define the //-bias to be A/x = /WB-oniy — A*s : Ms is the 
distance modulus from the E-fit method and ^t m s-oniy is 
the distance from the m^-only method. Fig. [T2l displays 
A/i versus redshift for each scatter model. For the FUN- 
MIX and G10 scatter models, the maximum A/x variation 
is only ~ 0.01 mag over the redshift range of each survey. 
For the Cll models A/x varies by several hundredths over 
the redshift range of each survey. This pi-bias is caused by 
the bias in /3 combined with Malmquist bias that results 
in bluer SN with increasing redshift. Figure [JJJ] shows 
the mean fitted color value versus redshift, and the cor- 
responding jU-bias is A/x ~ A/3 x (c fit ) where A/3 is the 
bias in /3 and (c fit ) is the mean SALT-ii color. The red- 
shift dependence on Aps is therefore directly proportional 
to the redshift dependence of (c ). Figure [TBI also shows 
the mean fitted stretch parameter (x\) versus redshift. 
However, since the corresponding bias on a is less than 
0.01, the resulting /z-bias is less than 0.01 x 0.3 ~ 0.003, 
and is thus much smaller than the bias from the color 
term. 

To summarize, these bias tests show that the G10 
model is internally consistent and that the wavelength 
dependence of the intrinsic scatter can be described with 
both methods. In the first method the scatter is de- 
scribed by the uncertainties in the SALT-ii light curve 
fitting model, and the m^-only scatter matrix is used in 
the SALT2mu stage. In the second E-fit method the light 
curve fitting model uncertainties arc set to zero and the 
full intrinsic scatter matrix is used in the SALT2mu fitting 
stage. The difference in results between these two meth- 
ods, as applied to simulations based on the G10 model, 
are negligible. However, if the simulated scatter model 
includes large anti-correlations such as those suggested 
by Cll, then these two methods are not consistent and 
the traditional m^-only method results in a significant 
redshift-dependent bias in the distance moduli (Fig. [P2|). 
More specifically, the bias is a result of incorrectly using 
the SALT-II (i.e., G10) uncertainties in the light curve fits 
of the samples generated with the Cll model. 

6.3. Estimating the w-Bias 



Nearby O SDSS-II 



SNLS3 




0.2 0.4 0.6 0.8 

redshift 



0.2 0.4 0.6 0.8 1 

redshift 




0.2 0.4 0.6 0.8 

redshift 



0.2 0.4 0.6 0.8 1 

redshift 



Fig. 12. — Average difference between the distance modufus 
from the mg-only fit and that from the E-fit, versus redshift. The 
simulated (true) scatter model is indicated in each panel. The sim- 
ulated samples are nearby (solid triangles), SDSS-II (open circles), 
and SNLS3 (solid circles). 
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Fig. 13. — Mean fitted SALT-II color vs. redshift for the Nearby 
(solid triangles), SDSS-II (open blue circles) and SNLS3 (solid red 
circles) simulation. Right panel shows mean fitted SALT-II stretch 
parameters (xi) vs. redshift. Although this simulation uses the 
G10 intrinsic scatter model, there is little difference for the other 
scatter model s. A lso note that the analogous plots in the bottom 
panels of Figs lU2l show a larger variation because of the more strict 
selection requirement of three passbands with S/N > 8. 



To estimate the bias on the dark energy parameter 
w, we fit the simulated Hubble diagram in a manner 
similar to that described in §8 of K09. The effect of 
Malmquist corrections is discussed later in £16.41 The w- 
bias is defined to be the difference between w obtained 
from the m^-only method and the E-fit method. Pri- 
ors are included from measurements of Baryon Acous- 
tic Osci llations from the SPS S Luminous Red Galaxy 
sample ([Eisenstein et al.ll2005t ) and from the Cosmic Mi- 
crowave Back ground temperature anisotropy measured 
from WMAP (jKomatsu et~all[200l . 

Since the simulated SN sample consists of 20 times the 
data statistics, the uncertainties on the priors have been 
reduced by a factor of a/20- We checked this procedure 
by splitting the simulated sample into 20 independent 
sub-samples and fitting each sub-sample with the nom- 
inal priors; the average of the 20 fitted w-valucs is in 
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TABLE 6 

a AND P FIT RESULTS FROM SIMULATIONS WITH DIFFERENT INTRINSIC SCATTER MODELS, AND UI-BIAS RESULTS a FROM COSMOLOGY FITS. 
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Intrinsic 














Malmquist- 


Scatter 


Scatter 






fitted 


fitted 




uncorrected 


corrected d 


Model 


Matrix 


xV^dof 




a h 


/9 C 


/3-bias 


Mi-bias 


Mi-bias 








SDSS-II 










FUN-MIX 


<T int 

m B 


= 0.087 


1.004 





114 ±0.003 


3.25 ± 0.04 


-0.01 ± 0.01 


±0.007 ± 0.006 


-0.022 ± 0.005 


G10 


int 
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±0.04 ± 0.01 


-0.002 ± 0.006 
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int 
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a Bias uncertainties are MC statistical uncertainties. 
b Simulated »sim =0.11. 
Simulated /3 S im = 3.20. 

d For m^-only method, correction is from simulation using fitted j3 value and G10 intrinsic scatter model. 
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good agreement with that from using the full simulation 
and priors with reduced uncertainties. The /3-bias and 
Malmquist-uncorrected iw-bias values are defined as the 
difference between the m^-only and E-fit methods, and 
they are shown in Table [5] The uncertainty on the w- 
bias is given by rms/v 7 ^), where rms is the dispersion 
among the 20 independent w-bias measurements. 

For the FUN-MIX and G10 models the w-bias is small 
(< 0.02) for all sample combinations. This small bias is 
expected since the fitted ft and distance moduli show a 
very small bias. For the Cll models, the w-bias is larger. 
Fitting the SDSS-II or SNLS3 sample alone, the w-bias 
is - 0.07 for the C11_0 model and - 0.12 for the Cll_l 
model. Combining the SDSS II + SNLS3 samples, the 
w-bias is reduced by a factor of ~ 2 for each of the Cll 
models. Including the nearby sample (nearby+SDSS- 
II+SNLS3) the w-bias is further reduced: 0.01 and 0.02 
for the C11_0 and Cll_l models, respectively. The re- 
duction in wj-bias as more samples are combined is a re- 
sult of the fortuitous circumstance that the mean SALT- 
II color for each sample is very close: —0.005, —0.015, 
and -0.005 for the nearby, SDSS-II, and SNLS3, re- 
spectively. To illustrate this point we have redone the 
analysis using a simulated nearby sample that has no 
Malmquist bias and a mean SALT-n color that is 0.06 mag 
redder than the data; the resulting ui-bias on the com- 
bined nearby+SDSS-II+SNLS3 sample increases to 0.05. 

6.4. Monte Carlo Correction for Malmquist Bias 

In £l6.3l the ui-bias is determined without accounting for 
differences in the Malmquist bias correction. While the 
true Malmquist bias should not depend on the analysis 
method (m^-only vs. S-method), here we show that the 
evaluated Malmquist bias, using the fitted value of ft, 
does indeed depend on the analysis method. In addition, 
using the evaluated Malmquist bias reduces the /i-bias 
shown in Fig. [T2] In the discussion below, "MqSIM" 
refers to the simulation used to evaluate the Malmquist 
bias correction. 

For the S-fit method the MqSIM uses the correct 
model of intrinsic scatter and the correct a and ft pa- 
rameters. In principle the intrinsic scatter matrix would 
have to be translated into a wavelength-dependent scat- 
ter model for the simulation, but we have not preformed 
this step. For the m^-only method the MqSIM for each 
sample uses the G10 model and the fitted (biased) value 
of ftm from Table [6] This procedure is used because it 
closely reflects the procedure in previous analyses. The 
distance-modulus corrections are from a second-order 
polynomial fit to /ifit — A*true versus redshift: /xtrue is the 
true distance modulus in the MqSIM, and iXfit is the dis- 
tance computed from the SALT2-fittcd parameters (color 
and stretch) and the simulated values of a and ft. 

Figure [14] shows the /i-bias versus redshift with the 
Malmquist correction applied. Overall the redshift- 
dependent bias is smaller than for the uncorrected dis- 
tances in Fig [12l but the bias is still significant for the 
Cll models. The bias is zero for the G10 model be- 
cause the fitted light curve model (G10) corresponds to 
the correct model of intrinsic scatter. The Malmquist- 
corrected w-bias results are shown in the last column of 
Table [5] Compared to the bias from the uncorrected 
Hubble diagrams, the iu-bias is typically smaller when 
the Malmquist correction is applied. 



While using an incorrect model of intrinsic scatter can 
lead to a significant bias in the Hubble diagram, this 
bias is somewhat reduced by simply applying a simu- 
lated Malmquist bias correction using the fitted value of 
ft in the simulation. The bias reduction depends on the 
intrinsic scatter model, and we cannot rule out increased 
sensitivity on other systematic effects. 



o 

f s 

I 
>• 



e 



o 

i o 

>- 

!= 

o 

i -o 

=L -0 

-0 
-0 
-0 



A Nearby O SDSS-II 



SNLS3 



FUN-MIX 



o- 



G10 



0.2 



0.4 0.6 

redshift 



1 o 0.2 



0.4 0.6 0.8 

redshift 



iO- 



C11 



■-o-o 



-o-o- 

* ^ 



C11 1 



0.2 0.4 0.6 O.f 

redshift 



1 o 



0.2 0.4 0.6 0.8 1 

redshift 



Fig. 14. — Same as Fig. 1121 but after applying Malmquist-bias 
correction. 



7. CONCLUSIONS 

We have used high quality SDSS-II and SNLS3 data 
and simulations to show that SN la intrinsic brightness 
variations include wavelength dependent variations re- 
sulting in a color dispersion of ~ 0.02 mag. A broad 
range of simulated intrinsic scatter models (Table [3]) arc 
roughly consistent with the following photometric ob- 
servables: Hubble scatter, dispersion in B — V — c', and 
photo- z residuals. These models include the G10 model 
that is dominated by a coherent term and has only posi- 
tive wavelength correlations, and the Cll model that has 
a small coherent term and large anti-correlations. 

We have used these intrinsic scatter models in high- 
statistics simulations to test the standard procedure of 
adding a constant distance- modulus uncertainty (<7i n t) 
to the measured uncertainties so that xV-^dof = 1 for 
cosmology fits to the SN la Hubblc-diagram. The con- 
stant CTi n t assumption is valid if the light curve fits in- 
clude model uncertainties with the correct wavelength 
dependence of the scatter. If these model uncertain- 
ties are not correct, such as using the SALT-ii uncertain- 
ties to fit a simulated sample with large anti-correlations 
in the scatter (Cll model), then significant biases can 
appear in the Hubble diagram. For the specific sim- 
ulation tests reported here the distance-modulus bias 
varies by up to 0.05 mag over the redshift range of 
the each survey, although this bias is roughly halved 
after applying Malmquist bias corrections. For the 
combined nearby+SDSS-II+SNLS3 simulated sample, 
w hich corresponds clos ely to the 472 SNe la analyzed 
in iSullivan et al.l (|2011[ ). the w-bias is only ~ 0.02, in 
part because the mean color for each sample is very 
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nearly the same. While this w-bias is well below the tota l 
systematic uncertainty reported in lSullivan et al.l ([201 1[ ) 
(dw sys t = 0.06), there is no assurance that future anal- 
yses on larger data sets will benefit from the fortuitous 
cancellations. For example, replacing the biased nearby 
sample with an unbiased sample increases the w-bias to 
0.05, comparable to the current systematic uncertainty. 

We have also shown that applying a simulated 
Malmquist bias correction, based on the fitted (biased) j3 
value, may reduce the bias from using an incorrect model 
of intrinsic scatter. However, we urge caution in relying 
on this apparent "free lunch" because we have not fully 
explored if this naive correction increases the sensitivity 
to other systematic effects. 

It is important not to interpret our results as proof that 
such biases exist in current SNIa-cosmology results; we 
have shown, for example, that this intrinsic scatter bias is 
negligible if the G10 model is correct. However, until the 
intrinsic scatter correlations can be better constrained, 
our results suggest that an additional systematic uncer- 
tainty should be included. 

Although the intrinsic scatter models considered here 
depend only on wavelength (except for the KRW09 mod- 
els), the true behavior could include a dependence on 
epoch, color, stretch, and redshift. The impact of the 
scatter models on the SALT-ii training is under investi- 
gation and will be published later. To obtain more ro- 
bust systematic constraints on cosmological parameters, 
we encourage additional studies to measure or constrain 
the nature of intrinsic scatter. 
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APPENDIX 

SIMULATION OF NEARBY SN IA SAMPLE 

Here we describe the simulation of the nearby (z < 0.1) SN la sample corresponding to the 123 nearby SNe la used in 
iConlev et al.l (|2011| ). While the nearby SN data ar e from several surve ys and filter sets, we simplify the simulation by 
considering only the CFA3-Kcplcrcam light curves (jHicken et al.ll2009f) ' using the UBVr filters. The i band is dropped 
because it is outside the valid wavelength range of the SALT-ii model. Since the CFA3-Keplcrcam simulation is a proxy 
for the entire nearby SN sample, we simulate the correct nearby-SN statistics corresponding to about half that of the 
SNLS3 sample. 

Since we do not have the observing conditions (mainly PSF and sky noise) needed for the SNANA simulation, we 
adopt a different strategy. The basic idea is to use each observed SN to define an observational sequence. The observed 
redshift and time of peak brightness are used along with the cadence. A random SALT-ii stretch and color are chosen 
for each SN, and the simulated S/N for each epoch is essentially scaled from the observed S/N. More technically, we 
artificially fixed the PSF to be 0.8" and then for each epoch compute the sky noise needed to simulate the observed 
uncertainty; this strategy allows generating nearby SNe in exactly the same manner as for the SDSS-II and SNLS3. 
With little knowledge of the spectroscopic selection criteria we cannot simulate the selection bias from first principles. 
We therefore modified the population parameters in Table Q] so that the simulated color and stretch distributions 
match those of the nearby sample. For the generated color distribution the peak-probability value depends on redshift, 
—0.02 — z/2, and the Gaussian width parameters are cr_ = 0.06 and cr + = 0.10. For the generated stretch distribution, 
the peak-probability value is reduced to 0.2. 

The intrinsic scatter matrix is determined for the nearby sample with the same procedure used on the SDSS-II 
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and SNLS3 simulations, and using SALT2mu we have checked that the fitted j3 values are in good agreement with the 
simulated input, /3sim = 3.2. We have also examined the analogous data/MC comparisons in Figs[T|[2l and find equally 
good agreement. 
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